Librerias

Primero incorporamos las librerias necesarias para poder ejecutar el presente notebook.

# Librerias
#install.packages('gganimate')
#install.packages("heatmaply")
#install.packages('gifski')
#install.packages('png')
#install.packages('units')
#install.packages('sf')
#install.packages('transformr')
require(tidyverse)
require(geosphere)
#<<<<<<< Updated upstream
require(rvest)
require(gganimate)
require(lubridate)
require(heatmaply)
require(transformr)
require(gifski)
require(png)
require(ggplot2)
#=======
require(rvest) # Cargamos el paquete
require(lubridate)
#>>>>>>> Stashed changes

Importar datos

Posteriormente importamos los datos que queremos analizar y los formateamos.

# Llamamos a los DataSets

datos_2021 = read.csv('../TP_LaboDeDatos/Data/sna_abril_2021_fixed_encoding.csv', encoding = 'UTF-8', sep = ',')

# Concateno todos los datos
datos2= read_csv2('../TP_LaboDeDatos/Data/202109-informe-ministerio.csv')

datos2$Fecha_completa = strptime(paste(datos2$Fecha, datos2$`Hora UTC`), format = "%d/%m/%Y %H:%M:%S")

En esa celda importamos los datos de los aeropuertos.

# Utilizamos la tabla que se encuentra en 'https://en.wikipedia.org/wiki/List_of_airports_in_Argentina'
# para acceder a las variables de ciudad, provincia y coordenadas de cada aeropuerto 

aeropuertos_wiki = read_html('https://en.wikipedia.org/wiki/List_of_airports_in_Argentina')
elemento_tabla   = html_element(aeropuertos_wiki,'.wikitable')
aeropuertos      = html_table(elemento_tabla)

Particularmente, necesitamos un procesamiento especial para obtener las coordenadas de cada aeropuerto en un formato que nos sea util.

# Corrijo la columna de coordenadas 

separar     = strsplit((aeropuertos$Coordinates), '/') # Divide a los strings en los lugares donde haya '/' 
coordenadas = sapply(separar, function(x) x[3])        # Me quedo solo con el 3er tipo de coordenada
coordenadas = gsub('[^0-9,.,-]','', coordenadas)           # Elimino los caracteres que no quiero utilizar

aeropuertos = aeropuertos %>% 
  mutate(lat = as.numeric(substr(coordenadas, 1, 9)), long = as.numeric(substr(coordenadas, 10, 18))) # Separo a mano latitud y longitud (revizar si esta todo en orden)

aeropuertos = filter(aeropuertos, nchar(ICAO)>1)
aeropuertos = aeropuertos[order(aeropuertos$ICAO),]

Procesar datos

En la siguiente celda lo que hacemos es a cada vuelo le agregamos los datos relativos a sus aeropuertos de origen y destino (si los tenemos). Además, calculamos con ellos la distancia recorrida por el vuelo.

for(i in 1:length(aeropuertos$ICAO)){
  inds = datos2$Aeropuerto==aeropuertos$IATA[i]
  datos2[inds,c('ciudad_origen','provincia_origen','lat_origen','long_origen')] = aeropuertos[i,c("City served","Province","lat","long")]
  
  inds = datos2$`Origen / Destino`==aeropuertos$IATA[i]
  datos2[inds,c('ciudad_destino','provincia_destino','lat_destino','long_destino')] = aeropuertos[i,c("City served","Province","lat","long")]
}
  for(i in 1:length(datos_2021$ana)){
  inds = datos2$Aeropuerto==datos_2021$ana[i]
  datos2[inds,c('ciudad_origen','provincia_origen','lat_origen','long_origen')] = datos_2021[i,c("nam","cpr","x","y")]
  
  
  inds = datos2$`Origen / Destino`==datos_2021$ana[i]
  datos2[inds,c('ciudad_destino','provincia_destino','lat_destino','long_destino')] = datos_2021[i,c("nam","cpr","x","y")]
  
  inds = datos2$Aeropuerto==datos_2021$iko[i]
  datos2[inds,c('ciudad_destino','provincia_destino','lat_destino','long_destino')] = datos_2021[i,c("nam","cpr","x","y")]
  
  inds = datos2$`Origen / Destino`==datos_2021$iko[i]
  datos2[inds,c('ciudad_destino','provincia_destino','lat_destino','long_destino')] = datos_2021[i,c("nam","cpr","x","y")]
}

datos2 = drop_na(datos2)
datos2$distancia = distHaversine(datos2[,c("long_origen","lat_origen")],datos2[,c("long_destino","lat_destino")])/1000

datos2 = datos2 %>% 
  filter(`Aerolinea Nombre` != 0)

Dividimos los eventos del dataset de vuelos en aterrizajes y despegues.

datos2$fecha_hora = strptime(paste(datos2$Fecha, datos2$`Hora UTC`), format = "%d/%m/%Y %H:%M:%S")
despegues = datos2[datos2$`Tipo de Movimiento` == 'Despegue',]
aterrizajes = datos2[datos2$`Tipo de Movimiento` == 'Aterrizaje',]

datos2

Ahora tratamos de unir los aterrizajes y despegues con el objetivo de poder identificar todo el recorrido de un vuelo. (particularmente la hora de salida y de llegada, que son los datos que no nos dan directamente en el dataset, y ademas la velocidad media teniendo el tiempo de viaje)

# dividimos en despegues y aterrizajes


matched = left_join(despegues, aterrizajes, by= c("Aeropuerto" = "Origen / Destino", "Origen / Destino" = "Aeropuerto", "Aerolinea Nombre" = "Aerolinea Nombre", "Aeronave" = "Aeronave")) %>% 
  mutate(tdif = as.numeric(fecha_hora.y - fecha_hora.x, units='hours')) %>%
  group_by(Aeropuerto, fecha_hora.x, Aeronave, `Aerolinea Nombre`) %>%
  filter(tdif > 0) %>% 
  filter(tdif < 5) %>% 
  filter(tdif == min(tdif))

matched = matched %>%
  mutate(vel_media = distancia.y/tdif)

Preguntas

Concentración del mercado

Las primeras preguntas que nos hemos planteado son relativas al grado de concentración que tienen los vuelos.

Por esto, la primera pregunta que nos interesa responder es ver cuales aeropuertos se conectan más y cuales menos. Además, notando que los casos donde el vuelo despega y aterriza en el mismo lugar son muy posiblemente vuelos de entrenamiento.

aers = unique(datos2$Aeropuerto)
mat = matrix(nrow = length(aers), ncol = length(aers))
colnames(mat) = aers
row.names(mat) = aers
for(i in 1:length(aers)){
  for(j in 1:length(aers)){
    mat[i,j]  =  sum((datos2$Aeropuerto==aers[i]) & (datos2$`Origen / Destino`==aers[j]))
  }
}
heatmaply(mat, Rowv = NA, Colv = NA,  col_dend_up=FALSE)%>% 
  layout(xaxis = list(side = "top"))

Ahora vamos a ver qué ocurre con el grafico anterior si lo pasamos a escala logaritmica (agregando 1 para evitar los -inf)

mat2 = log(mat+1)
heatmaply(mat2, Rowv = NA, Colv = NA,  col_dend_up=FALSE)%>% 
  layout(xaxis = list(side = "top"))

Como se puede ver, hay mayor cantidad de relaciones que se vuelven visibles y diferenciables gracias al cambio de escala.

Resalta el aeropuerto FDO como un importante centro de entrenamiento. Por otro lado, se observa como la conexión más fuerte es claramente entre Aeroparque y Bariloche.

La siguiente pregunta será relativa a la concentración entre distintas aerolineas.

aerolineas = unique(matched$`Aerolinea Nombre`)
apps = matrix(nrow = 1, ncol = length(aerolineas))
colnames(apps) = aerolineas
for(ar in aerolineas){
  apps[1,ar] = sum((matched$`Aerolinea Nombre`==ar))
}
apps = apps[1,order(apps, decreasing=T)]
#barplot(apps, las = 2, cex.names=0.4)
ejey <- list(title="Aerolineas", categoryorder = "array",
              categoryarray = rev(names(apps)), tickfont = list(size=4))
dfapps = data.frame(names(apps),apps)
dfapps = rename(dfapps,Aerolineas = names.apps., Cantidad = apps)
fig <- plot_ly(dfapps, x = ~Cantidad, y = ~Aerolineas, type = "bar") %>% 
 layout(yaxis = ejey)

fig

Y ahora probamos de repetir el ejercicio de expresar las distintas cantidades en escala logaritmica para ver el efecto.

appslog = log(apps+1)
#barplot(apps, las = 2, cex.names=0.4)
ejey <- list(title="Aerolineas", categoryorder = "array",
              categoryarray = rev(names(apps)), tickfont = list(size=4))
dfapps = data.frame(names(appslog),appslog)
dfapps = rename(dfapps,Aerolineas = names.appslog., Cantidad = appslog)
fig <- plot_ly(dfapps, x = ~Cantidad, y = ~Aerolineas, type = "bar") %>% 
 layout(yaxis = ejey)
fig

En este nuevo grafico se confirma la idea de que hay una enorme concentración de los vuelos en pocas empresas.

En este tercer grafico solo mostramos las más importantes para que sea entendible la información, además de que tiene una función más sintetica que los anteriores.

x = 8
apps2 = apps[1:x]
apps2[x+1] = sum(apps[x:length(apps)])
names(apps2)[x+1] = "OTRAS"
names(apps2) = paste(names(apps2),apps2)

fig <- plot_ly(data = as.data.frame(apps2), labels = names(apps2), values = apps2, type = 'pie', sort=FALSE)
fig <- fig %>% layout(title = 'Cantidad de vuelos por aerolinea',
         xaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
         yaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE)) 
fig

Tiempo en el aire

Las siguientes preguntas que nos hemos planteado fueron las relativas a cuando los aviones estaban volando y cuando no. Ya que nos parecia curioso ver cual era la cantidad de aviones que volaban simultaneamente, dicho de otra forma cuando se “cruzan”, independiente al destino que tengan

La primera fue decir exactamente en qué horas del año estuvieron más y menos aviones en el aire. Para eso primero hay que agregar una hora universal a los vuelos

inicio = as.Date("2021-01-01")
matched$despegueUniversal = as.integer(-difftime(inicio,matched$Fecha_completa.x, units="hours"))
matched$llegadaUniversal = as.integer(-difftime(inicio,matched$Fecha_completa.y, units="hours")+0.9999)
fin = max(matched$llegadaUniversal)+1
histo = matrix(nrow=1,ncol=fin)
histo[,c(1:fin)] = 0
for(i in 1:length(matched$despegueUniversal)){
  inds = matched$despegueUniversal[i]:matched$llegadaUniversal[i]
  histo[1,inds] = histo[1,inds]+1
  #print(histo[1,inds])
  #for(j in matched$despegueUniversal[i]:matched$llegadaUniversal[i]){
  #  histo[1,j] = histo[1,j]+1
  #}
}
#plot(1:length(histo[1,]),histo[1,], type="s")
ts= paste((inicio + 1:fin/24), as.character(1:fin%%24))
eje <- list(title="Hora", categoryorder = "array",
              categoryarray = as.character(ts), tickfont = list(size=5))
temp = data.frame(histo[1,],Fecha =eje$categoryarray)
temp = rename(temp,Cantidad = histo.1...)
fig <- plot_ly(temp, y = ~Cantidad, x = ~Fecha, type = "bar") %>% 
 layout(xaxis = eje)

fig

Para evitar el ruido y las pequeñas variaciones diarias y semanales, vamos a repetir el grafico suavizado por día y por semana.

#plot(1:length(histo[1,]),histo[1,], type="s")

histo2 = matrix(nrow = 1, ncol = length(histo[1,])-23)
length(histo2)
## [1] 6533
for(i in 24:(fin-23)){
  #print((i-23):i)
  histo2[1,i-23] = sum(histo[1,(i-23):i])/24
}
ts= paste((inicio + 24:fin/24), as.character(24:fin%%24))
eje <- list(title="Hora", categoryorder = "array",
              categoryarray = as.character(ts), tickfont = list(size=5))

temp = data.frame(histo2[1,],Fecha =eje$categoryarray)
temp = rename(temp,Cantidad = histo2.1...)
fig <- plot_ly(temp, y = ~Cantidad, x = ~Fecha,  type = 'scatter', mode = 'lines') %>% 
 layout(title = "Análisis con dia movil",xaxis = eje)

fig
#plot(1:length(histo[1,]),histo[1,], type="s")
lim = 24*7
histo3 = matrix(nrow = 1, ncol = length(histo[1,])-lim+1)
for(i in lim:(fin-lim+1)){
  #print((i-23):i)
  histo3[1,i-lim+1] = sum(histo[1,(i-lim+1):i])/lim
}
ts= paste((inicio + lim:fin/24), as.character(lim:fin%%24))
eje <- list(title="Hora", categoryorder = "array",
              categoryarray = as.character(ts), tickfont = list(size=5))

temp = data.frame(histo3[1,],Fecha =eje$categoryarray)
temp = rename(temp,Cantidad = histo3.1...)
fig <- plot_ly(temp, y = ~Cantidad, x = ~Fecha,  type = 'scatter', mode = 'lines') %>% 
 layout(title = "Análisis con semana movil",xaxis = eje)

fig

En el grafico con el promedio semanal movil se puede apreciar especialmente la caida en la cantidad de aviones en el aire durante la mitad del año.

Otra incognita que nos surgió fue saber cual es el horario en el que los pasajeros viajan mas, es decir, cual es el horario que registra en promedio mas cantidad de pasajeros para todos los aviones en el rango de esa hora.

require(gganimate)
horario_vuelo = datos2 %>%
  filter(datos2$Pasajeros > 0) %>% 
  group_by(horario = hour(`Hora UTC`)) %>% 
  summarise(pasajeros = mean(Pasajeros))

plot = horario_vuelo %>% 
  ggplot(aes(x = horario, y = pasajeros, color = pasajeros)) + geom_line() +
    geom_point() +
  transition_reveal(horario)+labs(title = 'Cantidad promedio de pasajeros a lo largo del dia', x = 'Hora', y = 'Cantidad promedio de pasajeros')+ shadow_mark()+ scale_x_continuous(breaks = round(seq(min(horario_vuelo$horario), max(horario_vuelo$horario), by = 1),1))
animate(plot, renderer = gifski_renderer(loop = T))

A la hora de graficar esta informacion, surgieron dos alternativas que visualmente se veian bastante informativas, y que nos parecio un detalle bueno para dejar ambos graficos a modo de mostrar lo que se charló en las clases acerca de la dicotomia a la hora de como comunicar los datos y qué se aprecia de cada grafico distinto.

plot <- horario_vuelo %>% 
  ggplot(aes(x = horario, y = pasajeros, fill = pasajeros)) + geom_bar(stat='identity')+transition_states(horario,100,5)+shadow_mark()+labs(title = 'Cantidad promedio de pasajeros a lo largo del dia', x = 'Hora', y = 'Cantidad promedio de pasajeros')+ scale_x_continuous(breaks = round(seq(min(horario_vuelo$horario), max(horario_vuelo$horario), by = 1),1))
animate(plot, renderer = gifski_renderer(loop = T))

Por ultimo, otra de las preguntas que nos nació de ver el dataset es cuales son los aviones mas usados por aerolinea, esto se reflejaria en la cantidad de observaciones de cada modelo. Para evitar datos repetidos, utilizamos el dataset

Al haber tantas aerolineas, se pierde mucha informacion, entonces nos pareció mas certero reducir el analisis centrado en las 3 aerolineas de las que hay mas datos, es decir Aerolineas Argentinas, Flybondi y Jetsmart Airlines y ver que clase de informacion obtenemos.

## me quedo con las 3 mas importantes

aerolineas_data= matched %>%
  group_by(`Aerolinea Nombre`) %>% 
  count(Aeronave) %>% 
  filter(Aeronave != 0) %>% 
  top_n(n=3, wt = n) %>%
  filter(`Aerolinea Nombre` == "AEROLINEAS ARGENTINAS SA" || `Aerolinea Nombre` == "FB LÍNEAS AÉREAS - FLYBONDI" ||`Aerolinea Nombre` =="JETSMART AIRLINES S.A." ) %>% 
  arrange(`Aerolinea Nombre`) 
  
aerolineas_data = aerolineas_data[order(aerolineas_data$n),]
aerolineas_data$Aeronave = paste(rev(as.character(1:length(aerolineas_data$Aeronave))),aerolineas_data$Aeronave)
plotaero = ggplot(data = aerolineas_data, mapping = aes(x = Aeronave, y = n)) +
  geom_bar(stat = 'identity', color = 'blue', fill = 'blue')+
  labs(y = "Cantidad de vuelos por aeronave") +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 3),
        axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        plot.margin=unit(c(1,1,1.5,1.2),"cm"),
        strip.text = element_text(face = "bold", color = "black", hjust = 0, size = 8),
        strip.background = element_rect(fill = "gray90", color = "gray")) + 
  facet_wrap(~`Aerolinea Nombre`, scales = "free_x")+scale_color_viridis()
  ggplotly(plotaero)

Por lo visto, los aviones mas usados son el Embraer 190, y los Boeing 737-700 / Boeing 737-800 en el caso de Aerolineas Argentinas, en Flybondi el avion mas usado tambien es el Boeing 737 en sus distintas versiones, y por ultimo en Jetsmart Airlines su unico modelo de avion es el Airbus A320, lo cual coincide con lo buscado en los respectivos sitios de cada aerolinea.